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** A computer  code,  based  on  a finite  element  algorithm,  was 
developed  to  solve  the  governing  equations  for  a three-dimensional 
flow  around  a spherically  blunted  cone  at  an  angle  of  attack.  The 
solution  is  valid  for  a boundary  layer  region,  including  the  flow 
separation  in  the  cross-flow  direction  in  the  leeward  area.  The 
required  initial  and  boundary  conditions  for  this  solution  are 
obtained  from  a finite  difference  computer  program  for  laminar 
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compressible  three-dimensional  boundary  layer  flow.  Extension  of 
the  analysis  to  include  the  turbulent  transport  using  the  mixing 
length  theory  and  a two-equation  turbulence  model  is  described. 

The  results  demonstrate  the  capability  of  the  finite  element  metho 
to  predict  the  three-dimensional  boundary  region  flow  with  cross- 
flow  separation  in  the  leeward  area  and  to  model  the  effects  of 
turbulent  flow. 
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ABSTRACT 


A computer  code,  based  on  a finite  element  algorithm , was 
developed  to  solve  the  governing  equations  for  a three-dimensional 
flow  around  a spherically  blunted  cone  at  an  angle  of  attack. 

The  solution  is  valid  for  a boundary  layer  region,  including  the 
flow  separation  in  the  cross-flow  direction  in  the  leeward  area. 
The  required  initial  and  boundary  conditions  for  this  solution 
are  obtained  from  a finite  difference  computer  program  for  laminar 
compressible  three-dimensional  boundary  layer  flow.  Extension 
of  the  analysis  to  include  the  turbulent  transport  using  the 
mixing  length  theory  and  a two-equation  turbulence  model  is 
described.  The  results  demonstrate  the  capability  of  the  finite 
element  method  to  predict  the  three-dimensional  boundary  region 
flow  with  cross-flow  separation  in  the  leeward  area. 
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1.0  SUMMARY 


Numerical  determination  of  the  complete  flow  field  description 
around  external  stores  configured  as  blunted  axisymmetric  bodies 
at  finite  angle  of  attack  is  presented.  The  finite  element  com- 
puter program  is  employed  for  solution  of  the  three-dimensional, 
compressible  parabolic  Navier-Stokes  equations,  which  describe  the 
complete  flow  field  including  flow  separation  in  the  cross-flow 
plane  in  the  vicinity  of  the  leeward  body  generator.  The  appro- 
priately required  initial  and  boundary  conditions  for  this  solution 
are  obtained  from  a finite  difference  computer  program  which  solves 
laminar  compressible  three-dimensional  boundary  layer  flow  around 
spherically  blunted  cones  at  an  angle  of  attack.  Extension  of  the 
analysis  to  include  the  turbulent  transport  phenomena  is  presented. 
Three  numerical  experiments  were  performed  to  establish  the 
requirements  of  finite  element  network,  finite  element  sizes,  step 
size,  and  computer  times.  The  results,  which  do  not  correspond 
to  any  specific  experimental  test  condition  or  geometry,  show  that 
the  numerical  solution  algorithm  does  compute  the  developing  flow 
field  on  a blunted  axisymmetric  body  in  a convergent  and  stable 


manner 


2.0  INTRODUCTION 


2.1  Background 

Solution  of  the  compressible  three-dimensional  boundary  layer 

equations  for  flows  over  arbitrarily  shaped  external  and  internal 

configurations  is  a current  design  and  analysis  requirement  for 

numerous  important  systems  applications.  It  is  only  in  the  last 

few  years  that  efficient  numerical  techniques  coupled  with  large 

high  speed  computers  have  become  available  thereby  making  solution 

of  these  problems  tractable,  e.g.,  see  the  Proceedings  of  the 

Conference  on  Aerodynamic  Analyses  Requiring  Advanced  Computers, 

19751.  Three-dimensional  boundary  layer  solutions  have  been  used 

to  aid  in  understanding  and/or  designing  both  internal  and  external 

flow  environments.  For  example,  three-dimensional  boundary  layer 

solutions  have  been  used  to  characterize  mixing  in  scramjet  com- 

2 3 

bustors  and  chemical  laser  optical  cavities,  Zelazny  et  al  ' 

Cebeci  et  al1  have  developed  a computer  program  which  is  optimized 
for  wing  geometries  whereas  Kendall  et  al1  have  developed  a program 
for  the  external  flow  field  about  general  shaped  configurations 

4 

for  which  real  gas  effects  are  important.  Blottner  has  discussed 
the  status  of  the  numerical  computation  of  3D  boundary  layers 
wherein  it  was  noted  that  one  of  the  major  problems  is  the  deter- 
mination of  the  inviscid  flow  conditions.  Several  solutions  have 
been  reported  for  blunt  body  shapes  (e.g.,  sphere-cone)  where 
curve  fit  relations  have  been  used  to  represent  the  pressure 
field  with  the  other  inviscid  flow  conditions  being  calculated, 
see,  for  example,  Vvedenskaya5  or  Popinski  and  Davis6. 

This  investigation  focuses  on  one  particular  aspect  of  three- 
dimensional  flows,  the  laminar  and  turbulent  boundary  layer 
developing  in  the  leeward  region  for  a sphere-cone  body,  Figure 
2-1,  where  cross  flow  recirculation  effects  are  important.  This 
region  is  characterized  by  a strong  viscous-inviscid  interaction 
and  as  shown  in  Figure  2-1  can  result  in  generation  of  separated 
flows.  The  solution  is  obtained  using  the  finite  element  method. 
This  method  has  been  applied  recently  to  problems  in  fluid 
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Figure  2-1.  Separation  Bubble  Flowfield  Model 
for  Tliree-Dimensional  Separation,  ref.  18 


Figure  2-3.  Circumferential  Pressure  Distributions, 
30%  Bluntness,  a = 10°,  ref.  18 


Figure  2*2.  Longitudinal  Pressure  Distributions, 
30%  Bluntness,  a * 10°,  ref.  18 


Figure  2-4.  Circumferential  Pressure  Distributions, 
30%  Bluntness, o = 10°,  ref.  18 
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mechanics  ' ' , and  a comparison  of  the  finite  element  and  the 

more  established  finite  difference  method  has  been  made  recently 

7 

by  Popinski  et  al  . 

In  solving  the  three-dimensional  laminar  boundary  layer 

equations  for  a body  at  angle  of  attack,  difficulties  in  obtaining 

solutions  at  the  leeward  symmetry  plane  have  been  encountered  by 

8 

many  workers.  Moore  showed  that  the  cross-flow  momentum  equation 
has  a unique  asymptotic  solution  at  certain  conditions  and  non- 
uniqueness occurs  otherwise.  This  indeterminacy  was  attributed 
by  Moore  to  a lack  of  previous  history  in  the  description  of  the 
fluid  entering  the  leeward  region  from  the  windward  surface  of 
the  cone.  Solutions  in  the  leeward  region  have  been  studied  by 
a number  of  investigators,  e.g.,  Moore®;  Cheng1®;  Vvedenskaya5; 


Libby11'  12 ; and  Dwyer13, 


Murdock  obtained  solutions  in  the 


leeward  region  which  were  not  known  previously,  and  determined 
conditions  for  which  the  boundary  layer  in  this  region  was 
independent  of  the  out-of -plane  flow.  He  investigated  cases  for 
which  the  complete  boundary  layer  solution  does  not  exist;  for 
this  case,  he  concluded  that  the  boundary  layer  model  has  a 
defect  which  results  in  discontinuous  derivatives  in  the  leeward 
plane. 

Similarity  three-dimensional  boundary  layer  solutions  obtained 
by  integration  around  the  cone  in  circumferential  direction  have 
been  reported  by  a number  of  investigators6'  15'  16 ' 17 . These 
boundary  layer  solutions  give  satisfactory  results  for  most  of 
the  body  surface;  however,  at  finite  angle  of  attack,  there  persist 
small  regions  centered  on  the  leeward  cone  generator  where  no  solu- 
tion is  possible.  This  region  of  flow  around  a cone  at  finite 
incidence  is  complex  and  sonic  speeds  in  the  cross-flow  direction 
can  be  obtained.  This  in  turn  can  yield  mixed-type  flow  field 
description  involving  different  types  of  equations  which  can  lead 
to  additional  difficulty.  Secondary  flow  separation  can  eventually 
occur  in  the  cross-flow  plane,  and  the  presence  of  the  reversed 
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flow  velocities  contributes  to  the  instability  of  initial  value 
computations.  Additionally,  the  leeside  boundary  region  is 
dominated  by  the  viscous-inviscid  interactions,  and  the  assumptions 
of  constant  transverse  pressure  and  conical  outer  flow  may  not  be 
valid. 

Thus,  this  leeside  region  enclosing  the  cross-flow  separation 

exhibits  distinctly  non-similar  behavior;  in  particular,  lateral 

derivatives  may  be  of  the  same  order  as  the  transverse  gradients. 

Therefore,  classical  boundary  layer  theory  fails  in  this  region, 

due  at  least  to  increased  boundary  layer  growth  and  strong  lateral 

gradients.  In  view  of  the  unique  features  of  the  separated  flow 

region,  the  governing  differential  equations  should  include  the 

lateral  diffusion  terms.  In  addition,  non-similar  external  free- 

stream  conditions  (from  theory  or  experimental  flow  measurements) 

showing  distinct  adverse  azimuthal  pressure  gradients  should  be 

used.  Detailed  features  of  this  type  of  flow  have  been  hypothesized, 

and  substantiated  experimentally,  for  both  sharp  and  blunted 

18 

circular  cones  by  Stetson  . Figures  2-2  to  2-4  show  experimentally 
18 

measured  surface  pressure  distributions  in  both  the  longitudinal 
and  circumferential  directions  for  the  blunted  cone.  Adverse 
pressure  gradients  in  the  leeward  area  are  well  defined,  and  points 
on  the  locus  of  separating  streamline  are  indicated. 

1 8 

These  experimental  measurements  indicate  that  the  separated 
flow  region  is  not  in  direct  communication  with  downstream  in- 
fluences (base  region)  and  that  a predominant  flow  direction  is 
uniformly  discernible.  A steady,  three-dimensional  equation  system 
is  derivable  from  the  full  Navier-Stokes  system,  under  these 
assumptions,  which  contains  three-dimensional  convection  and 
pressure  distributions,  but  discards  as  second  order  the  diffusion 
operator  in  the  predominant  flow  direction  (only) . This  parabolic 
equation  system  can  be  solved  using  downstream  marching  procedures 
associated  with  the  classical  boundary  layer  equations,  but  is  a 
two-dimensional  elliptic  boundary  value  problem  in  the  cross-flow 
plane.  The  most  prominent  additional  complexity  inherent  with 


this  "parabolic  Navier-Stokes"  system  is  that  the  static  pressure 
distribution  has  become  an  unknown,  and  methods  for  its  solution 
can  be  used  to  classify  various  investigator's  approach  to  their 
solution. 

Solutions  for  a modified  parabolic  Navier-Stokes  system  for 

the  sharp  cone  at  angle  of  attack  have  been  reported  by  Lin  and 

Rubin  . In  their  comprehensive  analysis,  the  equation  system  is 

solved  after  simplifying  the  transverse  momentum  equation  by 

balancing  the  transverse  pressure  gradient  and  centrifugal  forces. 

The  required  surface  pressure  distributions  are  obtained  from 

experimental  measurements.  Regions  of  separated  flow  in  the 

cross-flow  plane  were  numerically  predicted  in  the  leeside  of 

the  cone,  using  these  pressure  distributions,  for  angles  of  attack 

up  to  twice  the  cone  half-angle.  The  resolution  of  the  cross- 

flow  separation  in  the  leeward  region  is  well  illustrated  by 

examining  the  stream  surfaces  in  a cross-plane  and  the  distribution 

of  the  boundary  layer  growth  as  computed  by  Lin  and  Rubin  and 

20 

shown  here  in  Figures  2-5  and  2-6.  Lubard  has  also  performed 
an  extensive  analysis  for  a sharp  cone,  using  a solution  of  the 
parabolic  Navier-Stokes  equations  coupled  with  the  Rankine- 

20 

Hugoniot  relations  to  jump  the  bowshock.  For  this  analysis  , 

the  transverse  momentum  equation  was  utilized  to  determine  surface 

pressure  from  the  shock  solution  and  good  agreement  was  obtained 

1 8 

with  the  data  of  Stetson  . Solution  on  the  leeward  side  includes 
the  separation  bubble  which  is  clearly  visible  in  the  cross-plane 
distribution  shown  in  Figure  2-7  (from  Ref.  20) . Similar  techniques 
using  the  same  system  of  equations  was  used  by  Lubard  for  calcu- 
lation of  the  flow  over  a blunted  cone  at  an  angle  of  attack. 

The  initial  conditions  at  the  sphere-cone  tangency  plane  were 

21 

provided  by  a separate  time  dependent  program 

2.2  Objectives  and  ADoroach 

The  objective  of  this  study  was  to  develop  and  employ  a 
numerical  analysis  capability  which  describes  the  flow  field  on 
the  leeward  side  of  an  external  store  configured  as  a blunted 
axisymmetric  body  at  an  angle  of  attack,  e.g..  Figure  2-1.  In 

-6- 
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Figure  2-5.  Projected  Stream  Surface,  = 7.95,  0=10°,  Tw/T_  = 0.41 . 


.,  Sonic  Line: 


(uJ  + wJ)l/2/a  = 1,  (a)  Rex  = 3.6  x 10s , a * 8°,  (b)  Rex  = 3.6  x 10s, « = 10°,  (c)  Rex  = 

2.1  x 10s.  a = 20°,  ref.  19 


Figure  2-6.  (a)  Boundary -Layer  Thickness,  M„  * 7.95,  TW/TQ  * 0.41,  Rex  = 3.6  x 10s 

Experimental  Data  (Tracy  1963);  Present  Theory,  (b)  Displacement 

Thickness.  ■ 7.95;  x * 1.17  in.'Re  ■ 1.1  x 106,  or»  20°  (Above): A;. 

6*.  a * 8°  (Below); A; Axisymmetric  6*  , ref.  19 
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addition,  this  analysis  capability  would  be  sufficiently  general 
such  that  the  boundary  layers  could  be  laminar,  transitional  or 
turbulent  and  cross  flow  recirculation  effects  were  important 
and  could  not  be  neglected. 

The  approach  taken  was  to  merge  and  modify, ' where  necessary, 
existing  analyses  (computer  codes)  into  a viable  solution  tech- 
nique for  the  external  stores  flow  field  computation.  Specifically, 

the  flow  over  the  forward  part  of  the  blunted  body  was  computed 

,22 

using  an  existing  computer  program,  Popinski  , which  solves  for 
the  laminar  3D  boundary  layers  for  a sphere-cone  at  an  angle  of 
attack  where  cross  flow  recirculation  is  not  important.  Continua- 
tion of  the  solution  is  transferred  to  a second  computer  code  for 
conditions  where  the  assumption  of  laminar  flow  is  invalid  or 
cross  flow  diffusional  effects  become  important.  This  second 
computer  code  was  developed  from  a modification  of  the  finite 

element  analysis  used  in  studies  of  confined  (ducted)  three-dimen- 

2 3 

sional  laminar  and  turbulent  flows  ' . This  code  requires  a 

prescribed  external  pressure  distribution  which  as  discussed  by 

4 

Blottner  is  one  of  the  major  problems  with  solving  the  three- 
dimensional  boundary  layer  equations,  i.e.,  prescribing  the  inviscid 
flow  conditions.  In  this  study,  initial  and  inviscid  conditions 
are  used  to  perform  numerical  experiments  which  allow  the  influences 
of  turbulence  and  cross  flow  diffusion  to  be  isolated.  Three 
distinct  cases  are  considered  which  are  defined  by  the  combination 
of  initial  and  inviscid  conditions  employed,  Table  2-1. 

The  first  case  tests  the  ability  of  the  analysis  to  predict 

the  flow  field  on  the  leeward  side  downstream  of  the  sphere-cone 

junction  where  cross  flow  diffusional  effects  become  important. 

The  inviscid  data  are  based  on  the  pressure  distribution  obtained 

23 

using  the  method  of  characteristics,  Rakich  . The  second  case 
examines  the  effect  of  perturbinq  the  inviscid  pressure  field 
parametrically  to  induce  separation  in  the  cross  flow  plane.  The 
third  case  examines  the  effect  of  considering  transition  to 
turbulent  flow  on  the  leeward  side  by  "computationally"  tripping 
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the  flow.  These  calculations  are  essential  in  providing  the 
necessary  familiarity  of  the  sensitivity  of  theoretical  predictions 
to  critical  input  parameters.  It  is  recommended  that  other 
conditions  and  data  be  examined  in  subsequent  efforts  to  test  the 
generality  of  the  model. 

The  governing  equations  used  to  solve  the  flow  field  downstream 
of  the  sphere-cone  junction  are  given  in  Section  3.1.  The  finite 
difference  analysis  used  to  generate  the  initial  conditions  at 
this  junction  have  been  reported  elsewhere6  and  will  not  be 
repeated  here.  However,  key  features  of  the  forebody  flow  field 
solution  are  presented  in  Section  3.2.  Details  of  the  inviscid 
pressure  field  and  method  of  handling  external  boundary  conditions 
are  discussed  in  Section  3.3.  Various  turbulence  models  which 
have  been  examined  and  the  specific  approach  used  is  described  in 
Section  3.4.  A number  of  well  defined  test  cases  were  used  in 
the  course  of  this  study  to  establish  the  numerical  accuracy  of 
the  model  and  its  sensitivity  to  discretization  and  step  size. 

This  effort  included  establishing  the  relative  importance  of  the 
convective,  diffusion,  and  pressure  source  terms  in  the  momentum 
and  energy  equations  for  laminar  and  turbulent  flows.  Section 
4.0  describes  the  results  of  these  test  cases  and  the  results  of 
the  analysis  of  Cases  I-III  of  Table  2-1.  Conclusions  are  given 
in  Section  5.0.  Appendix  A presents  the  results  of  an  order  of 
magnitude  analysis  which  reduces  the  Navier  Stokes  equation  to 
the  parabolic  form  in  the  main  flow  direction.  Details  on  the 
finite  element  solution  algorithm  have  been  reported  in  Refs. 

2,3  and  33-37.  Those  modifications  which  were  required  to  extend 
the  existing  program  to  consider  the  subject  problem  are  discussed 
in  Appendix  B,  Finite  Element  Formulation.  The  method  used  to 
solve  the  continuity  equation  is  given  in  Appendix  C. 
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3.0  ANALYSIS 


3.1  Governing  Equations 

The  governing  equations  in  curvilinear  coordinates  for 
steady  state  conditions  are  given  as  follows. 


Global  Continuity 


drq  <h2h3pu)  + dq  (hih3pv>  + dq  (hih2pw)  * 0 


(1) 


x^  Momentum  Equation 
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x2  Momentum  Equation 
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Entha lov  Equation 


/ u 3H  , v 3H  w_  3H_  ) _ 

\ E^-  3x  3x2  h3  3x3  / 


Equation  of  State 


p = pRT 


Diffusion  Coefficients 


ue  = w + tpe 


Vi  = v (T)  Laminar 


T = Intermittency  function 


= Eddy  viscosity  (Section  3.5) 


u ■ Pe 
Pr  Pr„ 


= Turbulent  Prandtl  number 


Enthalpy-Temperature 
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Tg (H»  u,  v,  w,  Cp) 
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As  indicated  in  Section  2.0,  this  study  was  directed  towards 
spherically  blunted  circular  cones  at  an  angle  of  attack.  For 
this  geometry,  the  metric  coefficients  h^,  h 2,  h3,  and  the  required 
derivatives  are  given  in  Table  3-1  for  a body  oriented  coordinate 
system  as  shown  in  Figure  3-1. 


Examining  the  viscous  stress  tensor  and  performing  an 
order  of  magnitude  analysis  (Appendix  A) , the  governing  equations 
result  in  the  Equations  (2),  (4),  and  (5)  simplifying  to  the 
following  nondimens ional  form. 
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Table  3 . 1 

Sphere-Cone  Geometry:  Definition  of  Metrics 
For  Coordinates  Downstream  of  Sphere-Cone  Junction* 
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♦Body  oriented  conical  coordinate  system 
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In  this  study,  the  variation  of  pressure  across  the  boundary 
layer  is  neglected*  and  Equation  (2)  reduces  to 

= 0 (12) 

The  method  of  solution  of  these  equations  are  given  in  Appendices 
B and  C. 

3.2  Solution  of  the  Forebody  Flow  Field 

In  this  study,  an  existing  computer  program  developed  by 
22 

Popinski  was  used  to  obtain  details  of  the  flow  field  in 
Regions  I and  II  of  Figure  3-1  and  provide  initial  conditions 
along  surface  line  AB.  The  stagnation  point  solution  is 
determined  by  solving  the  similarity  equations,  obtained  by  a 
suitable  transformation  which  eliminates  dependence  on  the 
streamwise  coordinate  and  reduces  the  governing  equations  at 
the  stagnation  point  to  a set  of  ordinary  differential  equations. 
With  the  known  stagnation  point  solution,  solution  along  the 
windward  plane  is  obtained.  The  boundary  layer  solution  for  the 
stagnation  region  in  the  wind  coordinate  system  is  continued 
over  the  nose  beyond  the  stagnation  point  in  the  body  coordinate 
system.  A change  from  the  wind  coordinate  system  to  the  body 
coordinate  system  takes  place  at  a constant  value  of  the  surface 
coordinate  AB  (Figure  3-1).  Since  the  constant  value  Of  the 
surface  coordinate  in  the  wind  system  does  not  coincide  with  the 
constant  value  of  the  surface  coordinate  in  the  body  system, 
interpolation  of  the  computed  parameters  is  required. 

Examples  of  the  forebody  solutions  obtained  using  this  code6 
are  shown  in  Figures  3-2  to  3-5.  Figures  3-2  and  3-3  show,  for 
a sharp  circular  cone  at  angle  of  attack,  computed  values  for  skin 
friction  and  heat  transfer  for  the  classical  case  of  constant 
entropy  and  for  the  case  with  entropy  swallowing.  The  solution 
for  a spherically  blunted  conical  body  at  angle  of  attack  for  skin 
friction  and  heat  transfer  for  Regions  I,  II,  and  III  are  shown 
in  Figures  3-4  and  3-5.  These  results6,  based  on  solution  in 
' 21 

*See  Lubbard  and  Rakich  for  an  alternate  approach 
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Crocco  variables,  compare  favorably  with  results  obtained  in 

25 

Levy-Lees  variables,  Vatsa  and  Davis 

з. 3  Initial  and  Boundary  Conditions 

The  sphere-cone  body  under  consideration  as  shown  in  Figure 
3-6  identifies  the  coordinates  x^,  x2»  x3  and  the  corresponding 
velocities  u,  v,  and  w,  respectively.  The  flow  computations  in 
this  program  are  started  downstream  from  the  sphere-cone  junction 
point  (point  CJ,  Figure  3-7a) . The  initial  values  of  the  depen- 
dent variables  needed  as  starting  conditions  were  obtained  using 
the  approach  described  in  the  previous  section.  The  computed 
values  are  stored  on  an  external  storage  device  (magnetic  tape) 
and  interpolation  performed  to  obtain  values  at  the  finite  element 
node  points,  e.g..  Figure  3-7b.  This  interpolation  is  required 
since  the  mesh  used  in  the  finite  difference  forebody  analysis 
generally  must  be  more  refined  than  the  finite  element  mesh. 
Typical  initial  profiles  for  u,  v,  and  T in  the  starting  surface 
are  shown  in  Figures  3-8  to  3-10  for  selected  meridian  planes  for 
the  conditions  listed  in  Table  3-2. 

The  inviscid  edge  conditions  for  the  dependent  variables  to 
be  used  in  the  computations  as  boundary  conditions  have  been 
obtained  in  a form  compatible  with  the  finite  element  step  sizes. 
The  numerical  data  was  obtained  by  Rakich  using  a three-dimen- 
sional method  of  characteristics.  Figure  3-11  shows  the  circum- 
ferential distribution  of  the  inviscid  pressure  p,  velocities 

и,  v,  and  the  boundary  layer  thickness  in  the  starting  surface. 
These  data  are  based  on  pressure  distribution  for  a sphere-cone 
body  which  was  generated  by  the  method  of  characteristic 
(Rakich23) . At  this  stage  the  circumferential  pressure  distrib- 
ution used  in  the  computations  does  not  have  adverse  gradient 

in  the  leeward  area.  However,  once  the  program  is  made  opera- 
tional, the  pressure  distribution  will  be  modified  accordingly 
to  produce  flow  separation  in  the  cross-flow  direction  in  the 
leeward  area. 


Figure  3-7a.  Sphere-Cone  Geometry  with  Location  of  Starting  Conditions 
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Table  3-2 

Definition  of  Problem 


Cone  Half  Angle,  4>  /2 

C 

Nose  Radius,  (ft) 
Angle  of  Attack,  a 
Mach  Number,  M 

oo 


Reynolds  Number, 


PooUcoRn 

p 


10° 

1.0 

8° 

8 

8*105 


Wall  Temperature 
Prandtl  Number,  constant 
Ratio  of  Specific  Heats 


w 


*6To 

0.72 


1.4 
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Circumferential  Direction 


3.4  Turbulence  Models 


3.4.1  Background 

In  laminar  flows  the  boundary  layer  equations  represent  a 

closed  set  (equations  equal  dependent  variables)  since  the 

viscosity,  mass  diffusivity,  and  thermal  conductivity  are  known 

molecular  quantities.  However,  when  transitional,  turbulent, 

and  relaminarizing  flows  become  important,  the  additional  Reynolds 

stresses  are  introduced,  which  must  be  modeled  empirically  to 

complete  closure  of  the  equation  system.  The  approaches  taken  in 

modeling  the  turbulent  flux  terms  can  be  classified  into  two 

categories.  Category  I,  the  more  common  approach,  involves 

specification  of  eddy  coefficient  models  (eddy  viscosity  and 

turbulent  Prandtl  and  Schmidt  numbers) . Category  II  can  be 

defined  as  any  method  which  introduces  at  least  one  additional 

equation  to  describe  a turbulent  flux  or  a parameter  related  to 

a turbulent  flux,  e.g.,  turbulence  kinetic  energy  (TKE) , eddy 

viscosity,  length  scale.  An  excellent  review  of  the  various 

26 

models  for  turbulence  is  contained  in  Launder  and  Spalding  and 

27 

more  recently  by  Murthy 

It  is  significant  that  methods  of  Category  II  are  becoming 

increasingly  popular.  This  interest  may  be  attributed  to  the 

success  these  methods  have  had  in  accurately  predicting:  (1) 

28 

incompressible  boundary  layer  flows  (Kline  et  al  ),  (2)  transi- 
29  30 

tion  (Donaldson  ) , and  (3)  free  shear  layers  . The  ability  to 
model  a wide  range  of  flow  conditions  with  a single  set  of 
empirical  constants  is  the  main  advantage  of  Category  II  methods. 
This  is  not  surprising  since  the  physics  of  turbulent  flow  in- 
cluding upstream  history  effects  is  more  closely  modeled  than 
methods  of  Category  I.  However,  methods  of  Category  II  have 
certain  restrictions,  including:  (1)  requiring  empirical  relations 
to  model  the  higher  order  correlations  appearing  in  the  equations 
(data  needed  to  model  these  correlations  are  available  only  for 
limited  flow  conditions) , (2)  requiring  initial  profiles  for  the 
turbulence  parameters  being  considered,  and  (3)  an  increase  in 
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computer  computation  time  due  to  the  introduction  of  additional 
equations.  It  is  expected  that  with  the  continuing  interest  and 
effort  presently  directed  toward  methods  of  Category  II  that 
these  drawbacks  will  be  overcome  and  that  sophisticated  modeling 
of  turbulence  in  compressible,  three-dimensional  flows  will  be 
feasible  soon.  At  this  time,  description  of  turbulence  using 
methods  of  Category  I appears  most  plausible  because  they  have 
been  demonstrated  to  yield  solutions  to  practical  problems. 

However,  the  generality  of  our  approach  will  permit  using  methods 
of  Category  II  to  model  turbulence  as  these  more  sophisticated 
techniques  become  proven. 

The  location  of  transition  from  laminar  to  turbulent  flow 

is  dependent  on  such  parameters  as  surface  roughness,  freestream 

turbulence,  Reynolds  number,  and  pressure  gradients  (see  Morkovin31) . 

Characterization  of  this  region  is  important,  for  in  some  cases 

it  becomes  larger  than  the  laminar  and  turbulent  regions  and  the 

peak  heating  often  occurs  within  the  transition  zone  (Masaki  and 
32 

Yakura  ) . At  present,  there  is  no  generally  successful  technique 
for  predicting  transition.  The  transition  regions  have  generally 
been  specified  directly  as  input  into  the  solution  or  by  using 
empirical  correlations  based  upon  data  covering  a range  of  flow 
conditions,  see  Harris1. 


3.4.2  Approach 


At  this  time,  it  appears  most  practical  to  model  turbulence 
using  an  eddy  coefficient  hypothesis.  Following  Harris1,  the 
following  scalar  invariant  turbulence  model  is  used 
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(13) 


(14) 


(15) 


(16) 

(17) 


Equation  (7e)  is  used 
The  constants  k^ , 

values  of  0.435,  .09, 
respectively. 


to  obtain  the  effective  Prandtl  number, 
k^,  k^,  kg,  kg,  and  Npr  were  assigned  the 

26.,  0.01685,  5.,  0.75,  and  0.95, 


A second  turbulence  model  was  also  incorporated  into  the 

equation  system  which  retains  the  concept  of  an  eddy  viscosity 

but  relates  it  to  the  turbulence  kinetic  energy,  k,  and  dissipa- 
2 

tion  d as  e = Ck  /d  where 


(u'2  + v-2  + w"2)/2 


(18) 


The  length  scale  is  expressed  in  terms  of  k and  the  dissipation 
rate  for  turbulence,  d,  as 


l 


(19) 


The  turbulence  kinetic  energy  equation  is  introduced  to  determine 
k 


t » the  effective  stress  tensor 
u ■ time  averaged  velocity  vector 
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k,eff 


I 

I 

I 
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ii 

I! 

I 


Pr 


veff 


= effective  Prandtl  number  of  turbulence 
(empirical  constant) 

= dissipation  rate  constant 

= u + e 


In  many  flows  with  a predominant  flow  direction,  it  is  often 
sufficient  to  determine  the  turbulent  length  scale  by  equating  it 
with  the  characteristic  width  of  the  mixing  region.  However,  in 
many  complex  flows,  it  is  much  more  difficult  to  determine  a 
length  scale  due  to  the  complexity  of  the  flow  field.  To  overcome 
this  difficulty,  an  additional  conservation  equation  is  introduced 
for  the  dissipation  rate  of  turbulence 


In  the  following  section,  we  present  the  results  of  our 
efforts  to  modify  the  COMOC  computer  code  to  consider  turbulence 
and  cross  flow  diffusional  effects  and  discuss  future  plans. 


4.0  RESULTS 


4.1  Check  Cases 

The  parabolic  Navier-Stokes  equation  system  solved  using 
the  finite  element  solution  algorithm  may  be  routinely  employed 
to  consider  two-dimensional  problems  over  flat  plates  as  a 
special  case.  This  feature  is  important  since  changes  to  large 
computer  codes  easily  introduce  errors  which  often  can  only  be 
detected  using  existing  well  defined  solutions.  In  addition, 
the  accuracy  and  convergence  characteristics  of  solutions  to 
flow  fields  which  are  less  complex  than  the  3D  flow  problem  of 
interest  herein  can  be  used  to  gain  insight  into  optimum  mesh 
size  requirements,  i.e.,  conditions  which  minimize  computer  time. 
With  this  point  in  mind,  a compressible  laminar  boundary  layer 
developing  on  a flat  plate  was  successfully  used  to  establish 

2 

the  accuracy  of  each  change  required  to  modify  an  existing  code 
to  apply  to  the  three-dimensional  flow  over  blunted  axisymmetric 
bodies . 

As  a second  means  of  testing  the  accuracy  of  the  modified 

COMOC  code,  the  forebody  solution  at  the  sphere  cone  junction 

was  analyzed  in  detail  to  determine  the  relative  magnitude  of 

each  term  in  the  governing  equations.  Table  4-1  compares  the 

boundary  layer  thicknesses,  velocities,  pressures,  temperatures, 

viscosities,  and  pressure  gradients  between  a two-dimensional 

38 

compressible  boundary  layer  with  a known  analytical  solution 
and  the  three-dimensional  boundary  layer  of  interest  in  this 
investigation.  The  boundary  layers  on  sphere  cone  are  considerably 
smaller  than  the  flat  plate  case  due  to  the  highly  favorable 
pressure  gradient  upstream  of  the  sphere  cone  junction.  Figure 
4-1  shows  how  the  rate  of  change  of  momentum  is  distributed  across 
the  boundary  layers  in  both  the  streamwise  and  cross  flow 
directions  at  the  plane  which  will  be  used  to  start  the  afterbody 
flow  field  solution.  The  values  were  obtained  by  using  a finite 
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•Similarity  solution  from  Ref.  38  for  M = 5,  8 = 0.5,11  6/v  =830 

••Forebody  solution  for  conditions  listed  in  Table  3-2  . 


difference  technique  to  compute  each  term  for  comparison  with  the 
corresponding  terms  in  the  finite  element  code.  These  results 
were  then  used  to  aid  in  program  checkout  and  gain  insight  into 
the  dominant  mechanisms  governing  the  boundary  layer  development. 
The  two-dimensional  boundary  layer  is  also  shown  for  comparison. 

The  primary  difference  between  the  character  of  the  2D  and  3D 
flows  is  that  for  the  sphere-cone  body  the  boundary  layer  is 
being  accelerated  at  a much  higher  rate  (more  than  a factor  of 
10)  in  both  the  and  x3  directions.  The  larger  accelerations 
associated  with  the  three-dimensional  boundary  layers  result  in 
a stiffer  system  of  differential  equations  than  the  corresponding 
2D  problem  and  consequently  longer  computer  solution  times  are 
required  as  shown  below. 

4.2  Demonstration  Cases 

Table  2-1  showed  three  cases  used  in  this  study  to  examine 
the  performance  of  the  finite  element  numerical  solution 
algorithm  in  characterizing  three-dimensional  boundary  layer 
flows.  Items  of  particular  interest  were  computer  solution 
time  requirements,  the  number  of  finite  elements  needed  to 
assure  convergence  and  required  step  size.  The  initial  values 
for  the  flow  region  being  considered  are  taken  downstream  of 
the  sphere-cone  junction  point  as  shown  in  Figure  3-7a.  The 
mesh  size,  number  of  steps,  distance  marched  downstream  and 
computer  time  requirements  are  given  in  Table  4-2  for  all  three 
cases.  Calculations  made  with  additional  nodes  show  no  significant 
effect  (<  4%  differences  of  x * 330)  on  computed  values  downstream. 
Using  less  nodes  in  either  the  y or  | direction  resulted  in 
unstable  solutions.  The  length  coordinate,  x,  was  selected  to 
reflect  the  number  of  initial  boundary  layer  thicknesses  (at  4 « 0) 
which  have  been  inarched  downstream  of-  the  starting  location. 

Hence,  x - 333  corresponds  to  a distance  of  x ■ 6.88  ft  since 
6(x  » x^,  4 ■ 0°)  * .003  ft  and  x^  * 5.88  ft. 


S I 


The  edge  values  of  pressure  velocities,  temperature,  and 
boundary  layer  thickness  at  x = 330  are  shown  in  Figure  4-2  for 
Case  I.  No  significant  differences  over  the  values  at  x = 0 
are  observed  (see  Figure  3-11)  although  the  boundary  layer 
thickness  has  increased  approximately  a factor  of  twenty  over 
the  initial  values.  The  ten  node  points  used  across  the 
boundary  layer  were  found  sufficient  provided  transverse  length 
scale  (y)  was  normalized  by  the  boundary  layer  thickness. 
Computational  experiments  were  made  to  determine  that  over  65 
nodes  would  be  required  to  provide  the  same  degree  of  accuracy 
had  the  solution  been  obtained  in  physical  coordinates.  The 
velocity  (U  and  W)  profiles  computed  for  Case  I are  shown  in 
Figure  4-3. 

The  second  case  considered  represented  a perturbation 
about  Case  I where  the  only  parameter  changed  was  the  external 
pressure  distribution  as  shown  in  Figure  4-4a.  The  pressure 
variation  at  x = 330  was  found  to  induce  a cross  flow  recircula- 
tion pattern  at  x - 115  as  shown  in  Figure  4-4b.  This  result 
clearly  demonstrates  the  ability  of  the  finite  element  solution 
algorithm  to  predict  a smooth  and  stable  development  of  the 
recirculation  pattern  induced  by  the  adverse  pressure  distribu- 
tion prescribed  on  the  leeward  side  of  the  sphere  cone  body. 

The  third  case  studied  considered  the  effect  of  replacing 
the  laminar  flow  viscosity  model  with  an  effective  viscosity 
composed  of  the  sum  of  the  laminar  and  turbulent  viscosities. 
Equation  (20) . Three  types  of  turbulence  models  were  considered: 
(1)  mixing  length  theory,  (2)  turbulence  kinetic  energy.  Equation 
(20),  and  (3)  the  two-equation  turbulence  model.  Equations  (20) 
and  (21).  The  finite  element  solution  algorithm  was  modified 
to  consider  each  of  the  turbulence  models,  however,  computer 
time  limitations*  allowed  for  only  the  turbulent  mixing  length 

*Note  that  computer  times  for  solutions  from  X = 0 to  x ■ 170  with 
Equation  (20)  were  estimated  to  be  2.6  hours,  whereas  using  both 
Equations  (20)  and  (21)  result  in  computer  times  of  approximately 
3.3  hours  (IBM  360/65). 
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Figure  4-2.  Edge  Conditions  at  x * 330,  x = (x  - x0)  / 6 (*  - 0,  x = xQ),  Case  I. 
Laminar  Flow,  No  Cross  Flow  Recirculation 


(All  Cases) 


model  to  be  exercised  over  a distance  downstream  sufficiently 
large  to  make  comparisons  with  Case  I.  The  cost  of  the  turbulent 
flow  computation  is  almost  twice  that  of  the  laminar  calculation 
due  to  the  faster  mixing  rate  and  consequently  increased  rate  of 
change  of  the  dependent  variables.  Table  4-  compares  the 
predicted  skin  friction  distribution  at  x < 170  for  <j>  = 0,  90, 
and  180°.  As  expected,  the  "numerical  tripping"  of  the  flow 
at  x = 0 results  in  an  increase  in  the  skin  friction  by  as  much 
as  100%  over  the  laminar  value  at  x * 170. 
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Table  4-3.  Cf/2  Distribution*  Cases  I and  III 


Laminar  (Case  I)  Turbulent  (Case  III) 


<)>  (degrees) 

0 90  180 


*A11  values  are  multiplied  by  ^10' 


" '{T::  * 


5.0  CONCLUSIONS 

The  finite  element  solution  algorithm  used  in  the  COMOC 
code  has  been  extended  to  consider  cross  flow  recirculation 
effects  on  blunted  axisymmetric  bodies  at  a finite  angle  of 
attack.  This  recirculation  is  induced  by  an  externally  applied 
"surface"  pressure  distribution  which  produced  separation  on 
the  leeward  side  of  the  body.  Three  numerical  experiments  were 
performed  to  establish  the  requirements  of  finite  element 
network,  finite  element  sizes,  step  size,  and  computer  times. 

The  results,  which  do  not  correspond  to  any  specific  experimental 
test  condition  or  geometry,  show  that  the  numerical  solution 
algorithm  does  compute  the  developing  flow  field  on  a blunted 
axisymmetric  body  in  a convergent  and  stable  manner. 
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Appendix  A 

Order  of  Magnitude  Analysis 


The  solution  to  the  blunted  cone  flow  at  angle  of  attack  is 

based  on  an  approximate  system  of  equations  which  is  obtained 

from  the  steady  Navier-Stokes  equations.  These  equations  are 

derived  by  introducing  the  assumptions  that  the  streamwise 

derivatives  are  small  as  compared  with  the  normal  and  circumferential 
3 3 3 

derivatives  3^—  and  that  all  velocities  are  of  the 

order  o(l).  Terms  of  the  order  o(l)  and  0(6)  are  retained  in  all 
equations . 

The  divergence  of  the  viscous  shear  stress  tensor  for 
laminar  flow  is  given  in  curvilinear  coordinates  by  Equations 
A-l  to  A- 3 before  dropping  lower  order  terms. 

(V*T>  = 


h^hj  [xq  (h2h3Tn)  + drj  (hih 


3T12)  + 3x. 


. _1 3hl  . 1 3hl 

t12  h][HJ  SxJ  t31  h^ij  3x3 

1 3h2  t 1 3h3 

' t22  Kjq  3^  " t33 
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As  a result  of  the  order  of  magnitude  analysis,  the 
following  expressions  for  the  viscous  stress  tensor  terms  were 
obtained  and  used  in  the  parabolic  Navier-Stokes  equations. 
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Appendix  B 

Finite  Element  Formulation  and  Solution 

The  solution  of  the  viscous  flow  within  the  cross-flow 
separation  region  as  formulated  by  the  derived  parabolic  Navier- 
Stokes  system  is  obtained  using  an  extension  of  an  existing 
finite  element  computer  program.  The  finite  element  code  has 
been  applied  within  the  last  few  years  to  a family  of  elliptic 
and  parabolic  initial  boundary  value  problems  in  fluid  mechanics 
(Refs.  2,  3,  33-37),  and  their  relative  merits  as  compared  with 

7 

the  finite  difference  methods  have  been  recently  evaluated  . 

The  COMOC  code  is  based  on  a finite  element  numerical  solution 
algorithm  for  general  systems  of  nonlinear  elliptic  boundary 
value  partial  differential  equations  that  are  also  initial-value, 
i.e.,  parabolic.  The  general  differential  equation,  written  on 
a generalized  dependent  variable  q within  domain  R,  for  which 
the  finite  element  algorithm  is  established,  is  of  the  form, 

(puq),x  * (kq, ±)  , ± - (PUjq)^  + f(qqfi)  B-l 

subject  to  constraints  on  the  closure  3R,  identified  with  unit 
normal  vector  n^,  of  the  form  B(q)  = 0,  or  explicitly: 

B(q)  “ a^*q  + a^kq,jn^  - a*3*  = 0 B-2 

For  the  subject  problem  class,  the  generalized  diffusion 
term  will  contain  the  viscosity  coefficients,  y and  X,  and  the 
a^  are  user-specified  constraints  enforcing  all  forms  of 
physically  acceptable  boundary  conditions.  Table  B-l.  Solution 
of  Equations  9-11  is  accomplished  by  using  the  finite 
element  solution  algorithm.  Briefly  (Refs.  33  to  36),  approximate 
the  dependent  variable,  q,  within  M disjoint  interior  subdomains 
Rja  of  R,  defined  by  (x2,  x3,  x^eR^  x [xj,»),  where  the  union  of 
Rjg  is  R,  by  a series  expansion  of  the  form: 
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TABLE  B-l 

GENERAL  BOUNDARY  CONDITION  STATEMENT 


Boundary  Conditions 

.(1) 

.(2) 

.(3) 

No  Slip  at  Wall 

1 

0 

0 

Slip  at  Wall 

+ 

0 

Mass  Injection 

0 

1 

♦ 

Adiabatic  Wall 

0 

1 

0 

Specified  Heat  Flux 

0 

1 

+ 

Temperature  Dependent  Flux 

+ 

1 

+ 

Symmetry  Condition 

0 

1 

0 

♦User  specified  as  non-zero  to  enforce  desired  condition  level. 
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B-3 


I 

I 


f 

I 

II 


li 


0 

I! 


{l(x2,x3) IMQIXj) lm 


In  Equation  B-3,  the  elements  of  the  column  matrix  { Q } are  the 

m 


(unknown)  values  of  q at  the  nodes  of  Rm,  and  the  elements  of 
{L}  are  known  polynomials  (Equation  B-7) . 


The  natural  coordinates  {L}  of  any  point  in  the  discretization 
plane  are  expressed  implicitly  in  terms  of  the  coordinates  of  the 
nodes  in  the  local  coordinate  system  (Table  B-2) . The  finite 
element  solution  algorithm  applicable  to  elliptic  partial 
differential  equations  which  is  of  the  form: 


1 W(j>*  dx  = f Wf.  <♦*,♦?.  ,x.)dx  + i Wf.(<)>*)d a B-4 

^ kk  jRm  J3R_ 


m 


where  <p  is  the  approximation  to  the  dependent  variable,  f^  are 
specified  functions  of  their  arguments  and  contain  no  second 
derivatives  of  <p  , and  the  weighting  functions  W are  identical 
to  the  approximation  function  for  <)>  . Within  a finite  element 


m 


m 


B-5 


W = {L} 


B-6 


and  in  two-dimensional  plane  of  discretization 


1 - — + (-l-i)  H- 

V 'v  • \w 


(L) 


x_ 

X. 


x2  y3 


B-7 


and 
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TABLE  B-2 

IMPLICIT  DEFINITION  OF  SIMPLEX  NATURAL 
COORDINATE  FUNCTIONS 


TABLE  B-3 

INTEGRALS  OF  NATURAL  COORDINATE  FUNCTION 
PRODUCTS  OVER  FINITE  ELEMENT  DOMAINS 


D - Determinant  of  coeffident  matrix  defining  the 
natural  coordinate  tyttam.  tea  Table  B-2. 


*»  , 


I 

(TABLE  BA 

STANDARD  FINITE  ELEMENT  MATRIX  FORMS  FOR  SIMPLEX 
FUNCTIONALS  IN  ONE-  AND  TWO-DIMENSIONAL  SPAC  E 


Matrix*1*  Mitrix  ... 

Nam*  Function  Matrix  Evaluation'* ' •'*' 


(1)  Matrix  names  ara  a t digit  coda  covering  dimensional 1 ty , nonlinearity,  degree  of 
differentiation  and  ipocle)  Matrix  properties,  as  [a,  b,  c,  d,  e,  f]  where: 
a • A,  I,  C for  spaces  of  one-,  two-,  and  threa-dlmens  1 ons , 
b • number  of  coordinate  functions  appearing  in  Integral  or  isatrlx, 
c.  d,  e ■ (0,1)  loolean  counters  Indicating  (no,  yet)  differentiation  of  each  function, 
e or  f ■ S,  A,  A for  matrix  syaaetrlc,  antisymmetric  or  general. 

(I)  Syaaetrlc  Matrices  are  written  In  upper  trlongular  fora. 

(S)  A"  • 1/2  (X2P2)(X3P3) , the  plane  area  of  the  triangular  finite  element. 

1222  • the  x*  prlao  coordinate  of  node  2, 

X321  • the  Xj  prime  coordinate  of  node  3. 

(4)  »■  • length  of  side  for  boundary  condition  (*11222). 
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B-8 


I 

f 

f 

\ 

r 

li 

(1 

f! 

II 

n 

ii 

ii 

i; 

\ 

i 

i 

i 

i 


A A 

where  i and  j are  unit  vectors  in  the  local  reference  system. 

Integration  of  polynomials  over  domains,  which  is  required 
to  obtain  the  finite  element  formulation  of  the  partial  dif- 
ferential equations  is  straightforward  using  the  integration 
formula  shown  in  Table  B-3.  A series  of  reoccurring  integral 
expressions  is  presented  and  identified  in  Table  B-4 . 

The  computations  are  related  to  the  physical  conical  system 
consisting  of  an  annular  conical  ring  normal  to  the  surface  of 
the  cone.  The  conical  coordinate  system  and  the  physical  annular 
surface  are  shown  in  Figure  3-7b.  This  figure  also  identifies 
the  mesh  system  and  finite  elements  over  which  the  expressions 
are  evaluated.  The  curvilinear  coordinates  are  characterized 
by  the  shape  factors  h^,  h2,  and  h3  which,  of  course,  are 
evaluated  in  the  physical  conical  system.  For  computational 
purposes  by  finite  element  method,  the  annular  conical  element 
is  reduced  to  a rectangular  global  system  with  x2  and  x3 
coordinates  in  the  discretization  plane  and  the  streamwise 
coordinate  x^  normal  to  it.  The  relation  between  the  local  and 
global  coordinate  system  is  shown  in  Figures  3-7b  and  B-l . Form 
the  weighted  residual  of  Equation  B-l  by  substitution  and 
integration  over  R^  after  premultiplication  by  {L},  to  yield: 

{L}N(q*)dT  + XI  (L}B (q*) do  = 0 B-9 

jRm 

After  (a)  formation  of  Equation  B-9  for  every  element  Rm,  (b) 
Boolean  assembly  of  these  M equations,  (c)  integration  by  parts 
of  the  generalized  elliptic  operator  term,  and  (d)  identification 
of  the  Lagrange  multiplier,  X,  with  the  Reynolds,  Prandtl,  or 
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Schmidt  numbers,  as  appropriate,  the  global  finite  element 
solution  algorithm  is  established.  The  resultant  equation  system 
is  first-order,  ordinary  differential  and  is  written  on  the  x^ 
behavior  of  the  elements  of  {Q},  the  global  discretized  dependent 
variable  vector.  It  can  be  rendered  in  the  standard  form: 


jo]'  = 


B-10 


and  solved  for  Q using  an  integration  algorithm. 
Considering  the  x^  momentum  equation: 


L(u)  = 


with  boundary  conditions 


= 0 B-ll 


B (u)  = aj^u  + a2ku,ini  - a3 


B-12 


The  method  of  weighted  residuals  is  applied 

{L}L(u)dT  + if  {lJb  (u)  do  = 0 

jRm  i9\ 


B-13 


With  the  aid  of  the  Green-Gauss  theorem  and  assuming  that  the 
viscous  streamwise  derivatives  are  small  compared  with  the 

viscous  normal  and  circumferential  derivatives  - << 

} ■)  and  for  all  velocities  of  the  order  o(l)  (u,  v,  w = o(l)) 

one  obtains  the  following  form  of  the  x^  momentum  equation  after 
applying  the  method  of  weighted  residuals: 
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(L)  p 


u 3u 


v 3u  + w 3u  uv  1 

h2  3x2  FJ  3x3  hih2  3x2 


+ 


uw  . v2  3h2  w2  3h3~j  . 

3x3  h1h2  3xx  h!h3  8x1J  T 


3 i M_  3u 
3x2  h2  3x2 


* FJ  ^ (L>  Sj  %] 


Approximating  the  dependent  variable  q by  a series  expansion  of 
the  form 


vv 

* {L(x2,x3) }T{Q(x1) >m 

B-15 

where  the  elements  of  the  column  matrix  {Q}  are 

m 

values  of  q at  the  nodes  of  1^,  and  the  elements 
known  polynomial  shape  functions  (Equation  B-7) , 

the  unknown 
of  {l}  are 
with 

* 

“m  ’ 

{L)T{u(x1))m 

B-16 

<pu>m  ■ 

{l)t{rhooL, 

m 

B-17 

the  first  term  of  the  x^-momentum  equation  becomes  after 
evaluating  the  integrals  (Table  B-3) s 

K"  | {L}pu|£dT  " K"  { {L}{L}T{RHOU}in{L}T{u}^.T- 
1 jRta  1 ro 

{6 *2/2)  {2,2,1}  (2,1,2} 
{2,6,2}  {1,2,2} 
{2,2,6}. 


1 


Am 

60 


{RHOU}{u} 

B-18 


1 • 


Evaluating  all  terms  of  the  x^  momentum  equation  in  a similar 
way,  a first-order  ordinary  differential  equation  written  for 
uix^)  is  obtained  in  the  standard  form: 


= f ( 


[uj,*i] 


B-19 


which  can  be  solved  for  u at  the  downstream  station  using  any 
integrating  algorithm.  The  remaining  momentum  equations  and  the 
energy  equation  are  reduced  similarly  to  finite  element  form. 

The  first  term  of  the  energy  equation,  the  initial  value 
term,  assumes  the  following  finite  element  form  upon  integration 


f U>  f^“  dT 

Rm  3 X 


i-  {pu>T  ~ 

hj^  60 


III  111 
111  HI 


B-20 


Having  reduced  to  the  finite  element  form  all  terms  of  the  energy 
equation  this  way,  the  downstream  value  of  the  dependent  variable, 
the  total  enthalpy  H is  evaluated  by  integration  scheme. 

-—This  entire  system  of  equations  including  the  continuity 

equation  has  been  integrated  into  a single  code.  Extensive 
computational  exercises  of  the  flow  field,  including  a number  of 
variations  of  parameters,  have  been  performed.  A comprehensive 
documentation  of  the  pertinent  computational  results  is  presented. 


Appendix  C 

Solution  of  Continuity  Equation 


The  continuity  equation  for  a general  steady  state 
curvilinear  case  has  the  form 


^ (h2h3pU)  + j|j  (hjhjPV)  + (h^jPW)  ■=  0 


C-l 


After  introducing  the  finite  difference  expression,  the 
derivative  in  the  normal  direction  is  given  by 


jzr~  (h.h3pv) 
dx2  1 J i+l 


(h1h3pv)n  - (h1h3pv)n_1 


n-1 


dx. 


i+1 


jl-  (h2h3pu)  + 5^  <hiVw> 


n 


A+l 


C-2 


The  desired  value  of  the  normal  velocity  v”+1  at  the  n station  in 
the  normal  direction,  and  at  the  downstream  station  (£+1)  is 


(h^h3pv) 


n 

l+l 


(h^h3pv) 


n-1 

£+1 


53^-  (h2h3pu) 


n 


dx 


n-1 


C-3 


J £+1 


The  sum  of  derivative  in  the  expression  [ j is  evaluated 
appropriately  and  the  normal  velocity  is  determined  by  the  finite 
difference  method  in  each  column  marching  from  the  wall  toward 

the  edge  of  the  boundary  layer. 
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